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Abstract 

The dynamics of the avalanche width in the evolution model is described using 
a random walk picture. In this approach the critical exponents for avalanche 
distribution, r, and avalanche average time, 7, are found to be the same as 
in the previous mean field approximation but SOC appear at XcriUcai = 2/3, 
which is very close to numerical value. A continuous time random walk is 
studied numerically as a possible way to reconstruct in simpler concepts the 
evolution model. 
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INTRODUCTION 



Self Organized Criticality (SOC) theory studies non-equilibrium spatially extended sys- 
tems which show scale invariance on a large spatial or temporal domain. The phenomena 
which have these characteristics appear in astrophysics, geophysics, biological evolution, 
stock market, etc. [T|], 0, P . The study of this type of problems has emerged through 
various models; we will concentrate to the Bak & Sneppen model (BS) proposed as 
explanatory for certain aspects of the biological evolution. Its mathematical simplicity at- 
tracted numerical studies 0, 0, 101 and mean- field analytic treatments |Q. 

The model treats a number of species interacting on an one dimensional chain 
(a simple pictures of the food chain). Each species has assigned a scalar parameter, called 
fitness, with values in (0, 1) interval as a measure for the adaptability of species to the 
ecosystem. The dynamics starts from the site with the smallest fitness: new random values 
from (0, 1) interval are independently attributed to this site and to its two neighbors with a 
uniform distribution and uncorrelated from previous values then the cycle is repeated for a 
large number of times. For this system we define a A-avalanche ( < A < 1) as the number 
of steps between two consecutive configuration with all the fitness values greater than A. 
Numerically it has been found that the the stationary state exhibits scaling laws in ^ cxo 
limit for the correlation function of the active sites and avalanches distribution at A ~ 0.66 
0. We mention that the model shows also critical behavior in the versions with more than 
one dimensions or with slightly modified dynamics 1^, [0. 

In this paper we propose a new approximated solution. Our treatment brings in the 
dynamics of the avalanche width in a mean field approach. In the BS model the the avalanche 
width is a random variable which shows memory effect; we removed this effect by updating 
all the sites within an avalanche (see section II for more details), in this way we we kept 
track of the spatial extension of the avalanches, thing which is not possible in the infinite 
range approximation [Q. Even the approximation seems too crude we found the value of 
^critical = 2/3, which is very close to the numerical result. 
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The article organized as follow: Section I presents the general master equation for the 
fitness distribution and the derivation of the mean field equation found with probabilistic 
arguments in Ref . . Section II introduces our analytical treatment based on a mapping in 
a random walk problem. In this approximation we compute exactly the value of XcriUcai and 
the critical exponents r and 7 as they were defined in Ref. [^]. Section III presents a way 
to improve the method introducing a continuous time random walk description, numerically 
we show the improvement of the critical exponent r keeping XcriUcai fixed. The details of 
calculation are given in the Appendix. 



I. THE MASTER EQUATION. 

The BS model is completely characterized by the probability P{xi,X2, ■ ■ ■ , x^; t) to find 
the system in the state {xi,X2, ■ ■ ■ ,xn) at the time t given the initial distribution at t = 
0. Because there are not memory effects, the evolution of the system is described by the 
following master equation 

P(x.. X.. . . . . * + 1) = E / d.[M_M,.P.S; X,. . . . . x;. x;,. .... X.) 

i 

xP{xi, . . . , x[_i, x[, x'^j^i . . . , Xtv; (1) 

where periodic boundary conditions were assumed and Pstii] Xi, . . . , a;^_i, x^, x[j^^ . . . , xj^) 
is the probability to have activity at site i if the system is in the configuration 
(xi . . . , x',i_i^x\, x[_^_^ . . . xj\f). For the original one dimensional BS model 

Pstii] Xi,X2,...,xn) = Y[ ^i^j - ^i)^ (2) 
where d{x) is the step function 

fl, ifa;>0; 

9ix) = 

[0, ifx<0. 

At stationarity, integrating in (1) over X2, ■ ■ ■ iX^ we get easily the following relation: 

P„,(l, Xi) + Pae(2, Xi) + Pac{N, ^l) ' ^ = (3) 



where 

xP(xi, ...,Xn) 

is the the probabihty to have activity in site i when in site j the fitness has the value 
Xj. If in equation (|^) we try a stationary self-consistent mean field solution of the form 
p{xi, . . . , xn) = p{xi)p{x2) ■ ■ -P^xn) after some algebra we get 

(l - ]^)^''(^) + ]^^(^) + 3x - 3 = (4) 

with Q{x) = p{x')dx' . Equation was previously obtained in Ref. 0], in the N ^ oo 
hmit one finds p{x) =3/2 x e {XcriUcaU 1) and p{x) = when x e (0, Acritica/), Acriticai = 
1/3 whereas numerically Xcnucai ~ 2/3 0]. The statistical independence between the sites in 
the mean field solution allows the reduction of the problem at an one dimensional random 
walk on the positive semi-axis where the state n represents the state of the system with n 
fitness values greater than A. The solution developed in Ref. 0] gives the same XcriUcai as 
predicted by eq.(P and the critical exponents r = 3/2, 7 = 1. 

II. THE NEW APPROACH. 

We remember here that the size of an A avalanche is the number of steps between two 
consecutive events with no fitness below the the value A, so, it is a quantity characterizing 
time intervals. For a system of size N, with free boundary conditions, we define the avalanche 
width at a given moment t as the number of sites between the most left species with the 
fitness less than A and the most right species with the fitness less than A, the species between 
these two sites can have any value of the fitness. This is a quantity which characterize the 
spatial structure of our system. The width shows memory effects for the originally proposed 
dynamics, in the spirit of the mean-field approximation we approximate the evolution of 
the avalanche width with the following dynamics: at every step the species between right 
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and left extrema are performing uncorrelated movements and we introduce in dynamics the 
right nearest neighbor of the right extremum species in order to mimic the evolution. The 
complete randomness makes the movements of the two extrema completely equivalent and 
for this reason we have chosen to move only in one direction. In origin we also accept the 
double step. 

With this change the avalanche width is a random variable without memory effects on 
a discreet set of states which now can be extended to the entire non-negative semi-axis 
with the state zero corresponding to the state with no species below A and the state n to 
a realization of BS model with with n sites between the most left and the most right sites 
with the value of their corresponding fitness less than A. The transition matrix of the model 

has the following form: 

/(1-A)2 2A(1-A) A^ 

(1-A)2 2A(1-A) A^ 
(1-A)3 3A(1-A)2 2A2(1-A) 



P 





A^ 









;i-A)4 4A(1-A)3 3A2(1-A)2 2A2(1 - A) A^ 



V 



(5) 



The formulae for the matrix elements are: 

Poo = (1 - A)^ poi = 2A(1 - A), po2 = A^ poj = 0, j > 2; 
p,o = (1 - A)^+\ j>l; 

Pji = (j + 1)A(1 - Xy, J > 1; 

Pj2^jX\l-Xy-\ j>l; 

Pj+i^i-i if / > 2 and j < I — 1', 
otherwise 

with the convention that Pij is the transition probability from state the i to the state j and 
i,j e {0, 1, 2, . . .}. The distribution probability of avalanches is the first return probability 
distribution for this random walk and it can be written as 



Pji 



(6) 



p(n + 2) = J2PoiPi?Pio + IIpo2P2?P 



(7) 



i=l 



1=1 



where pj""* is the i,j element of the ra-th power of the matrix p obtained from the matrix 



p removing the row and the column zero and it is describing the evolution of the random 
walk outside of the origin. The first (second) term in the r.h.s. of eq.(|^ represents the first 
return probability, after n steps, when the initial step is single, (double). For a site different 
from the origin the forward step can only be single, as matrix p shows. 

We modify the first two columns of the transition matrix p such that to have the same 
elements on the diagonals of the p matrix. Keeping the closure relation J2jPij = 1 we 
produce the following matrix: 

/ (1-A)2 2A(1-A) 
(1-A)2(1 + 2A) 2A2(1-A) A^ 
(1-A)3(1 + 3A) 3A2(1-A)2 2A2(1-A) 



p 






A2 









(1-A)^(1 + 4A) 4\\1-Xf 3A2(1-A)2 2A2(1 - A) A^ 



V 



(8) 



The asymptotic behavior of the first return time distribution is the same for the both 
random walks described by the matrices p and p'. In fact in eq. (|^) we can go a step further 
developing p^^^ with respect to site 1 the remaining matrix from the second column and 
second row is identical with the p' matrix obtained from p' in the same way as p from p. 

3 /CO 



(n) 
Pit 



n ( H Pl2PhkPkl ) ( (1 - hi)Pl2P'2P + ^lipTl 



ni + ...nj 1 = 1 ^k=2 

— n — n' 



(9) 



with n = n' + rii + . . . rij and 6ij the Kronecker symbol. The terms in eq. (^) represent 
multiple returns in the site 1 before the last step to the origin. In the n — oo limit the 

~ (n) 

terms with n' and all the rij bounded but one have the same asymptotic behavior as p\^ 
because they are generated by the same matrix; the other terms will decay exponentially, 
due to P12 factor, or as a power of the leading term when there are two ore more unbounded 
exponents. U i = 1 n' has to be bounded to avoid the exponential decay. 

In the Appendix we present the computation for the generating function of the avalanches 
distribution probability ( [A15| ). In terms of the generating function R{z) = J2t^'^P{t) the 
average time for the avalanche distribution is 
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dz z=i d^ dz 2=1 



^From ( [A15| ) we obtain that the average time of an avalanche can be written as 



t 



A- A 



critical 



(11) 



when \\ — Xcriticail ^ 1, with the critical exponent 7 = 1 and the critical value of A, Xcriticai = 
2/3. We can also compute the asymptotic behavior of the avalanches probability distribution 
from the general formula 



+ dz' 



(12) 

20=0 



z=o 2m Jr {z - zqY^^ 

where F is an integration contour in the complex plane around zq which do not circle other 
poles. For A = 2/3 we found 

p{t) ^ t~^, t cx), (13) 

with T = 3/2, and exponential decay for A 7^ 2/3. In the previous equations the critical 
exponent r = 3/2 and 7 have values as obtained in the mean field solution whereas 
^critical = 2/3 is in extremely good agreement with the critical value of A found in numerical 
experiments [0, 0, [0. 

In language of Markov chain one can say that A = 2/3 is the transition point between 
persistent states (A < 2/3) and transient states (A > 2/3) [§. Nevertheless A is not a 
dynamical parameter for BS model, it introduces an "observational window" for a certain 
variable which we may choose from the set of statistical variables compatible with the 
dynamics of the BS model. SOC appears when there is at least one statistical variable with 
events at all scale lengths. In our approach A-avalanches are bounded to origin for A < 2/3 
and they escape to oo for A > 2/3. At A = 2/3 we have the peculiar stationary state in 
which the average time of avalanches is diverging, therefore there are events on the all time 
scales. 
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III. IMPROVING MEAN FIELD 



A significant difference between random walk proposed in the previous section and the 
BS model consists in the fact that in the latter the system will spend a characteristic number 
of steps in a given state because the activity can appear between the most left and the most 
right sites where the fitness is less than A. One possible way to improve our approximation is 
to promote the previous random walk to a continuous time random walk with inhomogeneous 
waiting time distributions, each waiting time distribution allowing for the persistence of a 
given size avalanche. The general equation for such a process can be written 

oo „J 

P^k{t) = 5.fce-^'* + ^ / c,e-^'V.,Pifc(t - t')dt' (14) 

where Pik{t) is the probability density to have the walker in state k at epoch t if at t = it 
was in state i, Pij are the elements of the p' matrix (^) and is the characteristic waiting 
time in site i and it represents the average life time for an avalanche of size i in BS model. 
Intuitively the average time of an avalanche is a function of the average number of sites with 
fitness less than A which is increasing with the avalanche width; at criticality we propose a 
behavior c^^ ^ for i > and Cq^ = a with a a given constant. 

An avalanche is now defined as an off time interval from the origin, whose probability 
distribution is independent of a 0. The avalanche distribution function can be expressed 
as 

oo 

Pav{t) = Pi^)Pn{t) (15) 
n=l 

where p{n) is the probability of an n steps excursion out of the origin and it is the first return 
probability distribution for the random walk defined in section II; Pn{t)dt is the probability 
of the first return to the origin in the interval t,t + dt after n steps. Intuitively we may say 
that for the A < 2/3 the exponential decay of p{n) will prohibit the long time avalanches 
and the average off time will be finite 0, while at A = 2/3 there is a qualitative change; 
even Pn{t) is decaying exponentially the scale invariance of p{n) for large n leads to critical 
behavior. 
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We have performed numerical simulation for the continuous time random walk at crit- 
icality for four scaling law of the parameters Cj, q = Xi^ with x = 1, 1.5, 2, 2.5. The 
numerical values for the critical exponent r (Table |) decreases monotonically as x increases. 
This behavior is intuitively clear, the avalanches tends to last longer if the characteristics 
life times are growing faster and r ^ 1 if x oo. 



IV. CONCLUSIONS 

We have proposed a new approach to Bak & Sneppen evolution model based on the 
dynamics of the avalanche width, the critical exponents are equal to those found previously 
in the mean field solution Q, in fact they are universal properties of the one dimensional 
random walk, but the value XcriUcai = 2/3 is very close to the numerical results reported in 
Ref. 0], [0]- Thereby we believe that the dynamics of the avalanche width is carrying useful 
information on the critical behavior for this model. The structure of the generating function 
( |A15| ) is intimately connected with critical behavior, the branch line appearing in ^ cxo 



limit generates the algebraic decay for the probability distribution of the avalanches. A 
generating function with a finite number of poles will lead, through formula (|12D, to an 
exponential decay of the avalanches probability distribution. The continuous time random 
walk picture allows for a more carefully analysis of the avalanche structure and it improves 
the r critical exponent keeping XcnUcai at the same value; it also gives an intuitive decom- 
position of the the algebraic decay distribution of the avalanches in a convolution of Poisson 
distributed events. 
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APPENDIX: 



We present the detailed calculation for the generating function of the avalanches proba- 
bility distribution. For this propose we shall use the special form of the matrix fl obtained 
from the matrix -p' removing the line and the row with index zero. This matrix has equal 
diagonal elements and we can write it as a linear combination of one-diagonal matrices /j 
defined as follow: 

/ i > 0, 

{h)u = { ' " (Al) 
Sk,i+i i < 



/q being the identity matrix. /^From the definitions ( |A1| ) we can compote the commutator 
T(*) = - for z > we have 

(T(^))h = (A2) 

and for ij > we have the property 

lilj = h+j- (A3) 

The matrix p' can be expressed as: 

oo 

P' = A^Ji + 2\\1 - A)Jo + + 2)A'(1 - \y+'I-^ 

i=l 

oo 

= A'/i + 1)(1 - \yi-i = X^hA (A4) 

i=0 

where A = I]^o(^ + ~ Xyi-i- Using eq. ( |A3|) it is easy to compute the n-th power of 
this matrix 



j=0 



J 



^From eq. (^) we can compute the commutator T„ = — A„Ji which has only the first 



column non zero 



{Tn)u=i^''^l^^yi-\)%i (A6) 
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consequently, 



(4T„),i= P^'+Z ^Vl-Ay+". (A7) 



J +n 

All the previous mentioned properties lead us to the relation 



n-2 



ihAr = I^A^ - ^ /,T,(/iA)"-^-2. (A8) 

i=0 

Eq. (|A8|) imply the following equation for the generating matrix G{z) = J^iZoi^'^-^^^T^^ 

oo 

G{z) = F{z) - J2 liTA'^'-^'^z'+'Giz) (A9) 

i=l 

where F{z) = Y^^ihA^z^ ■, z complex number. The sums which are appearing in eq. ( |A9D 
can be done in the following way: 

k=i k=i\ ^ + ^ / 

oo -| 

oo 1 Jk+j oo 1 „ 3k+j-l 

where = (1 — A)A^z. We can perform the summation if \ri^/{r] — 01 < ^^is set is 
not empty for < < 2/3-y/l/3. There is a annulus with inner radius and external radius 
obtained from the the positive solutions of the equation [r + ^y — r = 0; more than that, one 
of the roots of the polynomial —f]^ + V ~ ^ is inside of the minimal integration contour F for 



< ^ < 2/3yl/3 and the other two are outside of the maximal integration contour. Using 
the above mentioned properties of the matrices {Ik} (|A3|) we can compute the elements of 
the matrix /^A" which appear in the expression of the generating matrix F{z) : 

(^"''")'. = l) (1 - ^)""'"'. (^-"4")- = (, *' l) (1 - ^)""'- 

{LA")y = (InA")ij-i 3 > 1- 

Eq.(|^ shows that we need to compute only the first two rows in the generating matrices 
G{z) and F{z). The general formula for these matrix elements of F{z) is 
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2n 



n=j-l 

The previous series can be summed following the same computational path as in eq.( |Aiq ). 
If j = 1 we get 



for j > 1 we have the following expression: 

F (n = ^ ^ / 

(1- A)J-i27rUr -r/3 + r/-e' 
F2j{z) = Fij^i{z), j > 1, because Jn^" has equal elements on diagonals, and by direct 



calculation 



In all the above formula the contour F is the same as the that one used in eq. (|AlO|) and 
= (1 — A)A^z. Solving eq.([A9|) we obtain for the first two rows the solutions in terms of 
previously computed functions ui{z),U2{z), Fij{z), F2j: 



l + uAz) 



Residue theorem allow us to compute the generating functions in term of the third solutions 
of the polynomial —ff + Tj — ^^ VsiOy ^^at solution which lies inside of the integration contour 
r in the above integrals. 

1 

GlM) = (1 _ ^y^, ^'^^V ^ > 1' (All) 



^^^■(0 = -(T^V^l^ - 2)!^ , > 1; (A12) 
with ^2 = A2(l - A)^ and 

^ 1-^V3 _ (l + zv/3)(-27^ + V729e-108)^/^ 

22/3(_27e + V729e- 108)1/3 6 2V3 ' ^ ^ 



12 



^From eq. (|^) for the avalanches probabihty distribution one can write the generating 
function: 

oo oo 
R{z) = (1 - Xfz + Z^poi Gu{z)pio + Z^Pq2 G2i{z)pio 

i=l i=l 

^ -Fu{z){l-\Y+\l + {i + l)\) 



'l-\fz + 2\{l-\)z'Y.T 

i=i ^ 



+ ui{z)' 



+z'\' f:{F,,{z) - l'^'] ^ Fuiz))il - \y^\l + + 1)A). 



(AM) 



The series which are appearing above can be summed and the closed expression for gener- 
ating function reads: 

- 2^^)r,3K)^(l + 2A + (1 + SXM(? + + A^f^^))- (A15) 

The study of asymptotic behavior of the avalanches probability distribution is easily obtained 
by studying the behavior of the derivative of Gij{z) and G2j{z) functions. The summation 
appearing in R{z) changes only the amplitude for the leading term. It is crucial for the 
algebraic behavior the existence of the branch line in the complex plane for the generating 
functions Gij{z) which we choose on the real axis from (108/729)A^(1 — A) to oo. 
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TABLES 



TABLE I. The numerical values of the critical exponent r for four values of x 



X r 

1 1.27 it 0.01 
1.5 1.19 ±0.01 

2 1.14 ±0.01 
2.5 1.01 ±0.01 
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